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Abstract 


This thesis addresses the optimal control of wing rock problem, which is described 
by nonlinear equations of motion. A procedure is utilized to obtain optimal con- 
trol law for a nonlinear system, in which a positive performance index is minimised. 
The optimality equations are obtained from Hamilton- Jacobi-Bellman equation, which 
minimises this performance index. The positive definite closed loop Lyapunov func- 
tion is assumed to have same matrix form as performance index. Earlier, this method 
had been used for obtaining the solution of second order system, but in this thesis, 
it had been extended to third and fifth order systems. Both optimal nonlinear and 
linear feedback cases have been evaluated. It has been observed that both nonlinear as 
well as linear feedback control are equally efficient in controlling the wing rock motion 
of the model considered in this thesis at three angles of attack, 21.5 deg, 22.5 deg 
and 25 deg. A first-order actuator actuator has been incorporated and the effects of 
sideslip have also been considered. It has been shown that wing rock motion can also 
be initiated by sideslip even if sideslip angle is very small in magnitude. So, while 
designing the control system, this factor should also be given equal importance, which 
has generally been neglected in existing literature. It has also been observed that the 
simplified assumption made by Monahemi and Krstic that roll angle is twice of sideslip 
angle may not be applicable. 
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Chapter 1 


INTRODUCTION 


1.1 What is Wing Rock? 


The term “wing rock” describes the motion of certain types of aircraft such as fighters, 
which operate at angles of attack, either near to stall onset, or when flow seperation 
effects are evident by aircraft buffeting. This type of motion has similarity to a lightly 
damped Dutch-Roll mode, except that the wing rock motion is a limit cycle oscillation. 
Compared to Dutch-Roll mode the amplitude of the wing rock motion is solely a result 
of aerodynamic nonlinearities, while the amplitude of the Dutch-Roll mode response 
is determined by the initial conditions. 

Wing rock can also be described as lateral-directional instability for airplanes 
flying at high angles of attack. It is an uncommanded roll-yaw oscillation dominated 
by roll motion oscillating with a constant amplitude. It may be initiated either with 
a sideslip or during a zero-sideslip flight with some flow' asymmetries over the aircraft 
flying at high angles of attack. Once the asymmetric flow starts, a roll-oscillation 
amplitude will keep building up if roll damping is negative. The transient amplitude 
of wing rock will grow' gradually over some oscillation cycles because of instability and 
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negligible dihedral effect at low roll angle. 


Although the roll damping is negative at small roll angles, it is positive at larger _ 
roll angles for a sustained wing rock. Both the effective dihedral effect and positive 
roll damping via aerodynamic nonlinearities at large roll angles will gradually reduce 
the roll rate. As these restoring moments become stronger, the aircraft will reach 
a threshold roll angle and finally switch the rolling direction. Thus, wing rock is 
constrained to a finite amplitude oscillation through nonlinear roll damping. 

It occurs not only to low-aspect-ratio configurations, such as the F-4, F-5, F-14, 
X-29A, Gnat, Harrier, HP 115, and some re-entry vechicles, but also to high-aspect- 
ratio configurations, such as the Vari-Eze and some general aviation airplanes. Based 
on the extensive examination of available data, it can be concluded, as it has been done 
in Ref. [9], that wing rock is triggered by flow asymmetries, developed by negative 
roll damping, and sustained by nonlinear aerodynamic roll damping. 

Wing rock oscillations provide major limitations on the low speed maneuvers of 
those aircrafts which have highly swept wing planforms. At certain flight conditions, 
such as short landing on carriers at high angles of attack, these oscillations can cause 
safety problems that must be taken into account w'hile designing the aircraft control 
mechanism. Such mechanism can be established on the basis of some suitable mathe- 
matical model that is able to describe this highly nonlinear aerodynamic phenomenon 
as accurately as possible. The aerodynamic coefficients for this model can be obtained 
by measuring the forces acting on the aircraft, frequencies and their phase relation. 
On the other hand, in such experiments, these measurements are rather complex in 
nature because of the fact that the dynamic quantities involved are time dependent. 

So, in order to get the required data, Nguyen et al. [1] have performed an experi- 
ment in which they have successfully recorded roll angle oscillations during self-induced 
w’ing rock of a fiat plate delta wing with 80 deg leading edge sw T eep.- From this, they 
have calculated the resulting rolling moments. In their experiment, the rolling axis of 



free to roll model was mounted under the wing planform, resulting in a partial coning 
motion in addition to the pure roll oscillations. 

In continuation with these experiments, Levin and Katz [2] have performed a' 
similar experiment using two different models of aspect ratio 0.71 and 1.00. But this 
time their models have been mounted on a sting-balance apparatus and have been 
free to oscillate about their geometric axes, resulting in a pure roll oscillation. In 
addition to the roll angle, they have also measured normal and side forces and their 
phase relation. Theyhave conducted the testing of their models in two regimes, the 
static testing regime and the dynamic testing regime. They covered a wide range of 
angle of attack, angle of sideslip and roll angles to obtain the static test results. Based 
on these data, nonlinearities in the rolling moment slope have been detected. In the 
dynamic experiment, the wings have been mounted on a free to roll balance and the 
normal and side forces along with rolling angles have been then measured. Vortex 
patterns in the unsteady flow have been traced by helium-bubble flow visualisation, 
which have provided dynamic information about the vortex wake motion. This vortex 
wake motion plays an important role during the wing rock cycle. 

In the above experiment, it has been observed that the narrower wring with aspect 
ratio 0.71 has exhibited self-induced roll oscillations while the wing with the aspect 
ratio of unity, which has higher damping in roll coefficient, does not exhibit wing rock 
motion. During the free-to-roll tests a loss in wring average lift, relative to the static 
lift, has been detected for the same angle of attack. Apart from this, the oscillation 
frequencies and their variations with wrind tunnel speed and angle of attack, along 
with side force and normal force amplitudes, have been also recorded. 

The two major questions in regard to this phenomenon are: 

(1) What is the mechanism leading to the self induced oscillations? 

(2) What is the reason for the reduction in the average normal force? 

The explanation, as given in Ref. [2], based on the observed data is that these oscilla- 
tions are sustained due to the momentarily reduced distance between the up- moving 



semi-wing and leading-edge vortex originating from the same side. Thus, during the 
upward motion of this particular wing section a higher vortex suction force will result 
in a dynamic rolling moment, which is larger than the static value. This effect also 
causes the phase lag in the side force. On the other hand, vortex breakdown, has a 
damping effect on this motion, since flow visualisations have showed that the vortex 
burst is delayed relatively to the occurance of maximum roll angle. Furthermore, since 
under the dynamic conditions the wing upper surface is exposed to a more forward 
vortex breakdown, the average lift becomes less than the static value. In addition, 
during most of the roll cycle, only one leading edge vortex is present above the wing, 
whereas at the zero roll condition two of these vortices are present. Therefore, the 
reduction in the magnitude of the vortex lift, under the dynamic conditions, is obvious. 

One significant difference between the two experiments has been that, whereas 
Ngyuen et al. [1] have measured no wing rock for their 80 deg delta wing below 
a = 27 deg, Levin and Katz [2] have already measured wing rock at a = 20 deg 
for the same leading edge sweep, where a is angle of attack. This early wing rock 
occurance is probably, as the authors [2] suggest, caused by the centrebody axis that 
has been used in their model. Also the sideslips of the two experiments were different. 

Based on the information from these experimental results, Ericsson [3] has given an 
explanation of the underlying fluid mechanics leading to slender wing rock. According 
to him, two types of separation induced discontinuities occur for slender delta wing. 
One is caused by the leading edge vortices, for somewhat lesser slender delta wing. For 
very slender delta wings another type of discontinuity occurs before vortex breakdown, 
which is caused by asymmetric leading edge vortices. Vortex asymmetry occurs before 
vortex breakdown only for very slender delta wings. 

Ericsson [3] has observed that wing rock starts before vortex breakdown, and as 
has been suggested earlier, wing rock is associated with a loss of time-average lift. It is 
due to the fact that the lift-off of one of the leading edge vortices will cause a loss of lift. 
Thus, he has concluded that wing rock is caused by the vortex asymmetry and not by 



the vortex breakdown. He has explained the wing rock phenomenon in the following 
manner. At an a — 9 (apex half angle) combination where vortex asymmetry occurs, 
the wing half with the lifted-off vortex loses lift and dips down, rotating around the 
roll axis. As a result of increasing roll angle <f > , the effective apex angle 6 is increased 
and the vortex attaches again. This produces a restoring rolling moment, the positive 
aerodynamic spring needed for the rigid body oscillation in roll. Due to the convective 
time lag effect Ref. [3], the w r ing is dynamically unstable in roll untill the amplitude 
has reached the limit cycle magnitude. At this magnitude the damping on both sides 
of the discontinuity balances the undamping. 

Thus, the discontinuity which is introduced by the vortex asymmetry has all the 
characteristics that is needed for the limit cycle oscillation in roll. These character- 
istics, as he has observed, are lacking in vortex breakdown. If for some reason the 
vortex burst becomes asymmetric, the resulting net loss of lift on one wing half will 
cause it to “dip dow r n” . This increases (j) and thereby 0, causing the wing half to pen- 
etrate further into the vortex burst region. No switch to a restoring moment occurs. 
The opposite wing half gets out of the vortex burst region, generating increased lift 
which adds to the statically destabilizing rolling moment. Thus vortex breakdown has 
aerodynamic characteristics leading to roll divergence which are completely opposite 
, to those needed to cause slender wing rock. He also observed that the oscillations in 
roll damped down to zero amplitude if the 80-deg delta wing was yawed to 7 = 10 
deg at a = 27 deg. This is due to the fact that windward wing half has 8 > 15 deg 
leaving it outside of the boundary for asymmetric vortex shedding, whereas leeward 
wing half with 8 < 5 deg remains inside the region for vortex asymmetry. Thus, nei- 

wing half crosses the boundary and the wing-rock-inducing discontinuity is never 
mmuntered. 

In addition to this, Ericsson et al. [ 4 ] have made an extensive study of the ex- 
^©rinental results for advanced aircraft operating at high angles of attack. This has 

quite helpful in understanding the phenomenon of wing rock further. They have 
<Mws®rved that for aircraft with straight or moderately swept wings, wing rock occurs 
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first at moderate angles of attack. This type of wing rock is generated by the negative 
damping-in-plunge of the wing sections caused by dynamic airfoil stall. If the aircraft 
has long, slender nose, which often is the case, forebody vortices will generate wing 
rock at higher angles of attack where the wings are fully stalled. This type of wing 
rock has been observed in experiments with generic aircraft model, Ref. [4]. The flow 
mechanism causing the wing rock in this case is the downwash on the wing generated 
by forebody vortices. The X-29A aircraft experiences this type of forebody induced 
wing rock. When the inner portion of the wing has a highly swept leading edge, as in 
the cases of the F-18 HARV and X-31 aircraft, the high a wing rock is generated by 
direct interaction between the forebody vortices and the leading edge vortices from the 
inboard wing panels. Thus, from the above observations and a review of the available 
experimental results, Ericsson came to conclusion that wing rock of advanced aircraft 
will be driven by the vortices generated by forebody flow separation at high angle of 
attack, a > 30 deg. The key flow phenomenon driving the forebody-induced wring 
rock is the coupling between forebody crossflow separation and vechicle motion. More 
precisely, the rolling moment generates moving wall effects, that control the cross flow' 
separation on the forebody and the associated vortices provides the feedback mecha- 
nism responsible for the observed wing rock. The limit cycle oscillation results w’hen 
this driving rolling moment generated by the forebody vortices is balanced by the roll 
damping generated by the stalling wing. 


1.2 Mathematical Modelling of Wing Rock Motion 


After getting a good grasp of the wing rock phenomenon, now we proceed to analyse 
briefly the various valuable contributions of different scholars made in this field. 

Ross [5] has made a combined experimental and theoretical study of the nonlinear 
lateral oscillation experienced by a slender-wing aircraft at high angles of attack. 
Flight tests of HP 115 research aircraft have confirmed the expectation that the Dutch- 
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Roll oscillation is undamped at high angles of attack, and has also showed that the limit 
cycle develops, with steady amplitude at bank angle of about 30 deg. He has measured 
the stability derivatives along with the responses that has been obtained from the RAE 
flight dynamics simulater. From the digital computation of the equations of motion, 
he has observed that a wing rock limit cycle can exist when the nonlinearity is a cubic 
behaviour in either yawing or rolling moment with respect to the side slip angle. He has 
used the method of slowly varying parameters in order to demonstrate the limit cycle. 
In this method, the cross coupled nonlinear differential equations are reorganised into 
a single, high order differential equation, and the influence of nonlinear aerodynamic 
coefficients are analysed by stability criteria. 

An analysis is presented by Schmidt [6] using control theory concepts to show that 
aerodynamic hysteresis of the form of relay action can lead to lateral directional limit 
cycle motion. His analysis is an application of control theory, where the cross-coupled 
second order differential equations are cast into state variable form to yield a larger set 
of first order differential equations. This system is convenient for linear system analy- 
sis and quite adaptable to hysteresis situations from relay type inputs. The advantage 
of this technique is that physical concepts or principles become understandable and 
a broader insight is acquired. From his analysis, Schmidt [6] has made the following 
observations. 

(1) Limit cycle period does not differ significantly from the Dutch-Roll period, and 
does not change with the magnitude of roll moment hysteresis. 

(2) The limit cycle amplitudes are directly related to the magnitude of the roll moment 
relay action. 

(3) The sign of roll hysteresis parameter during switching can be determined from the 
requirement that a limit cycle with a dissipative system requires extraction of energy 
from the free stream. 

(4) During a limit cycle, the first harmonic of the roll rate response is approximately 
in phase with the sideslip angle rate while it is approximately 90 deg out of phase 
during a Dutch-Roll oscillation. 

(5) In addition to the phase differences between the limit cycle and Dutch-Roll oscil- 



lations, the relative amplitude ratios of bank angle to sideslip angle are different, for 
example, the ratio during limit cycle being about 8.5 times greater than that occuring 
during a Dutch-Roll oscillations. 

From the results of the experiments of References [1] and [2], a numerical simu- 
lation of these wind tunnel experiments has been developed by Konstadinopoulos et 
al.[7]. They have used unsteady vortex-lattice method to obtain aerodynamic loads 
and the equations governing wing rock motion have been integrated by a predictor- 
corrector scheme. The solution provides complete information about the motion of 
the wing and the flowfield simultaneously, fully accounting for dynamic-aerodynamic 
interaction. Their simulation predicts that the symmetric configuration of the leading- 
edge vortex becomes unstable as the angle of attack increases. This causes a loss of roll 
damping at small angles of roll. Consequently, small disturbances which are introduced 
into the flowfield grow, causing wing rock to develop. Also, from their simulation an 
analytic expression for roll moment can be developed, which is a third order polyno- 
mial in roll angle and roll rate that fits the numerical data extremely well from the 
start of the motion to the development of limit cycles or decay depending on the angle 
of attack. The model simulates hysteresis, aerodynamic damping and leading edge 
separation. It is not limited by angle of attack, twist, camber, or planform as long as 
vortex bursting does not occur in the near vicinity of the wing. 

Their numerical simulation predicts onset angles as well as the periods and ampli- 
tudes of the limit cycles in close agreement with the two experiments. Moreover, the 
simulation shows that the symmetric arrangement of the leading-edge-vortex system 
at zero roll becomes unstable. This asymmetry causes a roll moment to develop. But 
as the angle of roll increases, the change in the apparent direction of the free stream 
causes the direction of the moment and eventually, the direction of rotation to change. 

The aforesaid expression for roll moment has been substituted in the equation 
of motion by Nayfeh et al. [8] Thus, they have obtained a five term, third order 
polynomial in roll angle and roll rate. They have used this polynomial to construct 
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phase planes that show the locations of equilibrium positions, limit cycles, and domain 
of initial conditions that lead to the different possible motions. In addition to this, 
the use of this polynomial in the equation of motion makes it possible to obtain 
assymptotic approximations to the expressions for limit-cycle amplitudes and periods. 

Hsu and Lan [9] have developed another nonlinear aerodynamic model to inves- 
tigate the main aerodynamic nonlinearities causing wing rock. They have solved the 
resulting nonlinear flight dynamics equations of both one and three degrees of freedom 
with the help of Beecham-Titchener asymptotic method for the limit cycle amplitude 
and wing rock frequency. The latter are calculated iteratively in combination with a 
nonlinear aerodynamic method. Theoretical results of their model show that to sus- 
tain a steady state wing rock, the total aerodynamic roll damping must be negative 
at small bank angles. On the other hand, at large bank angles it must be positive. 


1.3 Control of Wing Rock 


There are three possible approaches to supress or prevent wing rock: 

(1) The first approach involves attaining the wing rock free capability through a de- 
tailed aerodynamic reshaping of the basic airframe configuration. 

(2) The second approach introduces maneuver limiting by addopting an alpha limiter, 
however, this approach degrades maneuverability. 

(3) The third approach employs stability augmented systems (SAS) or automatic flight 
control systems. These have become the most effective methods for attaining strong 
resistance to wing rock without degrading maneuverability. 

Luo and Lan [10] have presented a theoretical analysis to determine the optimal 
control input for wing rock suppression through a Hamiltonian formulation. The op- 
timality equations are analysed through Beecham- Titchener’s averaging technique 
and are numerically integrated by a backward differentiation formulas method. They 
have shown that the weighting factors in the cost function is related explicitly to the 
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system output damping and frequency. They have also shown that an effective way 
to suppress wing rock is to control the roll rate. System sensitivity is investigated by 
determining variations in system output damping and frequency with aerodynamic 
model coefficients. The results show that higher sensitivity corresponds to lower sys- 
tem damping. 

Monahemi and Krstic [11] have applied the theory of adaptive control of feedback 
linearizable systems for designing the control of wing rock motion with the extension 
of technique to include tracking. The adaptive regulation and tracking is achieved 
using a recursive design procedure known as backstepping with tuning functions. The 
procedure is such that at each step the subsystem is stabilised with respect to a 
Lyapunov function by designing an appropriate tuning function. At the final step one 
develops the feedback control and the update law for the parameter estimates. The 
designed, control and parameter update laws can track a prespecified model system 
with remarkable accuracy reaching steady state within a transient stage of any desired 
time length. Thus, without any degradation the required maneuvers can be maintained 
simultaneously with the suppression of wing rock. Such a performance requires the 
availability of sufficient rolling moment because of the lateral control. The design 
proceedure accomodates aerodynamic parameter estimation and updates online while 
the control system is in operation. 

Singh et al. [12] have employed an adaptive and neural control of the wing rock 
motion. Neural networks have architecture that has massively parallel interconnec- 
tions of simple “neural” processors. These networks have ability to “learn” a mapping 
between an input space and an output space. In this approach, however, the control 
cancels the nonlinearity of the main equation (the roll moment equation) only in a 
very limited sense, i.e., there exists no additional dynamics and / or an actuater control 
surface equation. 

Adaptive control systems developed in Refs. [11] and [12] assume that the wing 
rock model does not have disturbance input and that the structure of nonlinear func- 
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tions, except for the neural control, is known. In the presence of the external dis- 
turbance inputs and unmodeled nonlinearities, adaption law's of internal type used in 
Refs. [11] and [12] can cause divergence of controller parameters. 

Araujo and Singh [13] have designed a control system for the control of w'ing rock 
motion of simple slender delta wings which is based on the theory of variable structure 
model reference adaptive control (VS-MRAC), using only input and output signals. 
Their wing rock model can include disturbance input and unstructured nonlinear 
aerodynamic functions. In Ref. [13] only roll angle measurement is required for 
controller synthesis. For the derivation of the control law, it is assumed that the 
aerodynamic parameters and the structure of the nonlinear functions in the model 
are unknown. Moreover, it is assumed that the disturbance input due to w'ind gust is 
present in the system. Ref. [13] have shown that, in the closed loop system, the roll 
angle tracks a given reference trajectory, and wing rock suppression is accomplished. 

Shue et al. [14] have presented a procedure for optimising a state feedback con- 
trol law for a nonlinear system with respect to a positive performance index. They 
have employed Hamilton- Jacobi-Bellman equation to derive the optimality equations 
wherein this performance index is minimised. The closed loop Lyapunov function is 
assumed to have the same matrix form of state variables as the performance index. 
The constants of these matrices are easily determined so as to guarantee their positive 
definitenesses. The optimal nonlinear system is asymptotically stable in the large, as 
both the closed loop Lyapunov function and performance index are positive definite. 
It is also showm that the wing rock model using nonlinear state feedback is asymptot- 
ically stable in the large. Both optimal linear and nonlinear state feedback cases are 
evaluated. Even though the major emphasis is placed on the second order system in 
their paper, the procedure is applicable to any system of finite order. It is this pro- 
cedure, which has been utilised in this thesis, to obtain the optimal feedback control 
for wing rock problem along with some modifications like considering an actuator and 
sideslip effects, which have been neglected in Ref. [14]. 



Chapter 2 


METHODOLOGY 


2.1 Modelling of the Wing Rock Motion 

Currently there are several theoretical models proposed to describe wing rocking. Each 
model has been variefied onh r for a particular configuration. In fact, none of the models 
can predict the amplitude and period of wing rock motion for all configurations. 

In this thesis, the models developed in Refs. [8] and [12] have been used.The 
equations of motion of the tw T o References are compared here in order to get an idea 
of the variables on which the coefficients, used in Ref. [8] depend. Rate of change of 
angular momentum gives the following equations: 

L = I x p (2.1) 

or 

p = r < 2 - 2 > 

where p is the roll rate, L is the rolling moment and I x is the roll moment of inertia. 
Also from Ref. [12] we have 

P = [Cip (a) jlsinoc + C lp (a, 0){pb/ 2V)] (2.3) 

i x 
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where q, S, b and V are the dynamic pressure, wing reference area, wing span and 
freestream air velocity respectively. The coefficient Cip is the rolling moment stability 
derivative because of sideslip angle 13 and the coefficient Ci p is the rolling moment 
derivative because of roll rate p. On comparing Eqs. (2.2) and (2.3) we get 

L = qSb[C lfi (a)psina + C, p ( a, /3)(pb/2V)] (2.4) 

Now, rolling moment L depends on sideslip angle (3 and roll rate p. Therefore L can 
be written as 

L = Lp(3 + L p p (2.5) 

where L@ and L p are rolling moment due to sideslip angle and roll rate respectively. 
On comparing Eqs. (2.4) and (2.5) we gel, 

L@ = qSbCip(a)sina (2-6) 

L p = ,SbC lp (a,P)(pb/2V) (2.7) 

Thus, we can write 

p = j-J + 4 ? ( 2 - 8 ) 

The equation of motion used in Ref. [8] is 

4> = CiQ- C 2 ij> (2.9) 

where <j> is the roll angle and Ci is the aerodynamic roll-moment coefficient. 

4> = P (2.10) 

C, = pScL 2 c /2I x (2.11) 

C 2 = p x L c /I x U c (2.12) 

where p is the density of the air, c the chord, L the characteristic length, p x the 
damping coefficient for the sting bearing and U c is the characteristic speed. The dot 
in Eq. (2.9) implies that the derivative is with respect to the non-dimensional time 


r = ( u c /L c )t 


(2.13) 



(2.14) 


On comparing Eqs. (2.8) and (2.9) we get 

c 2 = -^C lp (aJ)b/ 2 V 

*X 

and from Eqs. (2.12) and (2.13) we find that 

Px = —qSbCip(a, (3)/ 2 


Again on comparing Eqs. (2.8) and (2.9) we get 

CiQ = ~C l0 {a)0sina 

*X 

where, 



Thus, on substituting the expression for q and C\ in Eq. (2.16) we get 


pScL 2 pV 2 Sb . 

-Q = —^7 — Ci 0 (a)/3sina 


2I X 


2I X 


On simplification we get 


V 2 

Ci = —Ci 0 (a)3sina 

cb 

The expression for Cj used in Ref. [8] is 


Cl — 0,\(j) + 02^ + 0-3 + C4<^ 2 <^ + CLglfrffi' 


(2.15) 

(2.16) 

(2.17) 

(2.18) 

(2.19) 

( 2 . 20 ) 


The value of a* depend on the angle of attack, as it can be observed from Eq. 
(2.19). They are given in Table (1). Eq. (2.20) fits the calculated moment through 
the transients and into the regime where steady state motion exists. This form for 
Ci has evolved from a much more general beginning, which included a total 12 terms 
in a fifth order polynomial. Seven of these terms have been systematically eliminated 
because of their small contributions. 

Substitution of Eq. (2.20) into Eq. (2.9) leads to 

<j> + u) 2 (f> = P\(j) + bxffi 4- p2<!? <j> + b2(j>4? (2-21) 

where, 

w 2 = — Cjflj (2.22 J 


t 4 
A — 



Hi = C\Q2 — C2 

(2.23) 

b\ = C1G3 

(2.24) 

Hi — 

(2.25) 

62 = C1G5 

(2.26) 


The point where hi is zero corresponds to the onset of wing rock. 


Table 1. Coefficients in the analytical expression for the 
roll moment, at different angles of attack 


a 

ax 

&2 

a 3 

d^[ 

0'S 

15 

-0.01026 

-0.02117 

-0.14181 

0.99735 

-0.83478 

21.5 

-0.04207 

0.01456 

0.04714 

-0.18583 

0.24234 

22.5 

-0.04681 

0.01966 

0.05671 

-0.22691 

0.59065 

25 

-0.05686 

0.03254 

0.07334 

-0.3597 

1.4681 


Now’ we have reached at a stage where wre must seek some approriate procedure 
to solve the above mentioned equation of motion, so that w-e may be able to observe 
the type of motion which the aircraft is executing, and latter develop some control 
system to control the undesirable motion. 


2.2 Control of Wing Rock Motion 


Before developing some procedure to solve our equation of motion obtained in the 
previous section, let us have a brief understanding of the control systems. The basic 
control system problem may be described by the simple block diagram show r n in Fig. 
(3.1). The main objective of the system is to control the variable c in some prescribed 
manner by an acuating signal e through the elements of controlled process. 




Error 

detector 


Command or 

AX 

Error or 

Forward path 

Output or 

reference input r 


actuating signal e 

elements 


controlled variable c 


Feedback path 
dements 


Figure 2.1: Basic block diagram of a feedback control system 

Control systems may be classified into two types depending on whether the con- 
trolled variable affects the actual input to the control system or not, i.e., whether some 
kind of feedback is used or not. When feedback is not used, the system is called an 
open loop system, when used, it is called a closed loop system. 

In an open loop system an input known as reference input is applied directly to 
the controller, and an output known as controlled output is obtained. Here in an 
open loop system the input to the controller is in no way affected by the value of 
the controlled output. Also, this type of system do not adapt to the variations in 
environmental conditions or to external disturbances. 

What is missing in the open loop control system for more accurate control is a 
link or feedback path between the output and the input of the system. In order to 
obtain more accurate control, the controlled signal must be fed and compared with the 
command or reference input, and an actuating signal proportional to the difference 
of the output and the input must be sent through the system to correct the error. A 
system with a feedback path is called a closed loop system. 

The analysis and design of feedback control systems are usually carried out by 
using either linear or nonlinear techniques. Strictly speaking, linear systems do not 
exist in practice, since all physical systems are nonlinear to some extent. Therefore, 
linear feedback” control systems are idealized systems. W hen the magnitudes of the 
actuating signals are limited to a range in which system components are considered 
linear, these components may be represented by their linear models and counterparts. 

zr. 




But when the actuating signals are extended outside the range of linear operations, 
the feedback control systems become nonlinear. 

We can now proceed forward to solve our equation of motion, under these two 
headings. 

2.2.1 The Open Loop System 

The equation of. motion of the aircraft is 

4> + up <j> = + b\ (P 4- /j, 2 (jP <f> + b 2 <f><jP (2.27) 

The above second order equation can be written as a system of two first order equa- 
tions. 

a?i = x 2 (2.28) 

±\ = x 2 = -uPx i + nxx 2 + b x xl + fJi 2 x\x 2 + hxixl (2.29) 

where 

X\ = 4> (2.30) 

To find the equilibrium positions, we substitute 

x x = f 2 = 0 (2.31) 

From this, it follows that 

*i = 0. ,±u/Jbi (2-32) 

The nonzero solution exists only when b x is greater than zero. It is observed that &i 
and /q change sign at the same angle of attack, Ref. [8]. At the angle of attack when 
bi is zero, two of these equilibrium positions are at infinity. As the angle of attack 
increases, these two point move towards the origin. 



2.2.2 The Closed Loop System 


Flight at high angle of attack, the exhibition of limit cycle behaviour of the motion 
variables, and the creation of large amplitudes for the motion variables lead to severe 
requirements for the flight control system that cannot be met using conventional pro- 
cedures. Nonlinear properties are essential features of the aircraft dynamics in these 
regimes. Techniques that are needed must be 

(1) computationally feasible, 

(2) retain the nonlinearities inherent within the vechicle dynamics, 

(3) provide physical and mathematical understanding of the important flying qualities. 
Thus, keeping the above factors in mind, the same procedure which has been used in 
Ref. [15], has been utilised in this thesis for optimizing a state feedback control law 7 
for our nonlinear system with respect to a positive performance index. Here also, the 
Hamilton-Jacobi-Bellman equation is used to derive the optimality equations so that 
this performance index is minimised. 

Let us consider a vector-matrix linear differential equation 

x = Ax + Bu (2.33) 

where x is the process state, u is the control input, and A and B are known matrices. 
We now seek a procedure to minimise a specified performance index (or cost function) 
expressed as the integral of a quadratic form in the state x plus a second quadratic 
form in the control u; i.e., 

T* 

J = J t [x T Qx + u r Ru]dt (2-34) 

where Q and R are symmetric matrices. Let us talk about this performance index in 
brief, before w r e attempt to find the optimum solution. 

First, we note that minimization of J also minimizes p* J where p* is any positive 
constant. So the problem is not altered by multiplying J by any positive constant. 
Often the constant 1/2 is used in front of J to simplify expressions resulting in other 
developements. 
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Second, regarding the limits on the integral, the lower limit t is identified as the 
present time, and the upper limit T* is the terminal time, or final time. The time 
difference T* — t is the control interval, or “time to go” . In our case, the terminal time 
is infinite. Here, we are interested in the behaviour of the process “from now on”, 
including the steady state. 

Finally, consider the weighting matrices Q and R. These are called the state 
weighting matrix and the control weighting matrix, respectively. The question of 
concern to the control system designer is the selection of the weighting matrices Q 
and R. 

In the performance or cost function defined as above, two terms contribute to the 
integrated cost of control: the quadratic form x T Qx which represents a penalty on 
the deviation of the state x from the origin and the term u T Ru which represents the 
“cost of control” . This means that the desired state is the origin and not some other 
state. The weighting matrix Q specifies the importance of the various components of 
the state vector relative to each other. The term u T Ru in the performance index is 
included in an attempt to limit the magnitude of the control signal u. Unless a “cost” 
is imposed for the use of control, the design that emerges is liable to generate control 
signals that cannot be achieved by the actuator and the result will be that the control 
signal will saturate at the maximum signal that can be produced. 

The above type of problem where performance index is of a quadratic form and 
the system is linear one, is known as linear quadratic regulator (LQR) problem. The 
Hamilton-Jacobi-Bellman (HJB) equation is 

H = ^x r Qx + ^u T Ru + ^^x (2.35) 

where V(x) is the Lyapunov function and denotes the derivative with respect to 
elements of x. 

Earlier, the (HJB) equation has been applied in solving (LQR) problems. LQR is 
guaranteed to have an optimal and asymptotically stable solution because the closed 



loop Lyapunov function is positive definite or positive semi-definite and its time deriva- 
tive is negative. 

problem, 

(2.36) 

(2.37) 

(2.38) 

(2.39) 

In fact, no matter whether or not the linear system or nonlinear system uses 
nonlinear feedback control, the system’s closed loop Lyapunov function has to be 
positive definite. As with the LQR problem, the performance index plays a very 
important role in the problem of optimal control of nonlinear systems. When a positive 
definite Lyapunov function cannot be found, this purpose usually can be achieved 
through changing the cost function. 

In this thesis, a procedure is employed that produces a closed loop Lyapunov 
function and a performance index that are positive definite, and the HJB equation is 
also satisfied. This means that the system is asymptotically stable in the large. 

"Consider a time invarient nonlinear system 


For the extension of the LQR problem to the nonlinear optimal control 
the following conditions must be satisfied. 

Condition 1: 

V(x) > 0, V{0) = 0 

Condition 2: 

• r rn . f 1 r _ 1 7V. dV(x ) . 1 

mvn{H} = mm < -x r Qx + -u T Ru-l — — ^-x > = 0 
(2 2 ox J 

Condition 3: 

x t Qx > 0 

From condition 2 we can find that 

V(x) = = -i[x T Qx + u r Ru] < 0 

if condition 3 is satisfied. 


x = Ax + F(x) v Bu 


(2.40) 



where A, B, x, u have already been stated, and F(x) represents the higher order term 
with respect to the state variable x. Therefore, /(x) w'hich is an element of F(x), may 
be w r ritten in the following form: 


/(*) = /a(x) + / 3 (x) + • • • + / 2 „_i(x); n > 2 


(2.41) 


Note that the subscripts of /(x), the performance index, the closed loop Lyapunov 
function, and the Hamiltonian matrix represent the order of the state variable x. We 
can define a performance index 

1 f°° r 


where 


ff(x) = 


] r°° 

J = ^ J [?( x ) + u T Ru]dt 


Ql,l Ql,2 • • • Ql,n 
Q2,l Q2,2 • • • Q2,n 


J- -j 

T f 

X 


X 2 


x n 




■ 


X 


X 2 


x n 


(2.42) 


(2.43) 


Qn,l Qn,2 • • • Qn,n 

with q(x)>0, and Qy = Qy when Also, R is assumed to be symmetric, orthog- 
onal, and positive definite matrix. Therefore, 


q{x) = q 2 (x) -h q 3 (x) -h ■ ■ ■ + ? 2 n(x) 


(2.44) 


where 


q 2 (x) = x r Q 1)1 x 
q 3 (x) = x t Q 12 x 2 + (x 2 ) r Q 2jl x 
94 (x) = x t Qi i 3 x 3 + (x 2 ) r Q 2j2 x 2 + (x 3 ) r Q 3 ,iX 


92n(x) = (x n ) T Q n>n X n 

The optimal control problem of the system in Eq. (2.40) is to find the state feedback 
gain 

u = u(x) (2-45) 

such that the performance index is minimised. The HJB equation is applied to find 
the optimal solution of the system, 

H = ^q{x) + ^u r Ru +'^y^[Ax + F(x) + Bu] = 0 


(2.46) 



where 


ar(x) 

dx 


dV(x) OV(x) OV(x) 


For optimal control, 


dx, ' dx 2 ’•••’ dx n 

dH 


d u 


= 0 


Thus on solving the above equation, we obtain the optimal control as 


u(x) = -R _1 B t 


dF(x) 


dx 


(2.47) 

(2.48) 

(2.49) 


On substituting Eq. (2.49) into Eq. (2. 40) wo get 


/f =5 9W -5^ BR “ Bl ' 


3F(x) 


dx 


+ ^^[ Ax + F ( x )] = ° (2-50) 


It is assumed that the Lyapunov function has the same form as the cost function (J), 

if, • 

>T r 

Pl.l Pi, 2 ••• Pl,n 

P 2,1 P2,2 • • • P2,r 


“ 5 


Pn,l Pn,2 ••• P n,n 


X 

v 2 


(2.51) 


paid P(x) > 0, where Py = Pj j when i j. therefore, 


where 


V(x) = F 2 (x) + V 3 (x) + • • • + V 2n (x) 

V 2 (x) = |x T Pi t iX 
V 3 (x) = |x r P 1>2 x 2 + |(xT P 2>1 x 


(2.52) 


P2n(x) = |(x n ) T P n , n X n 

Note that the performance index and the Lyapunov function are symmetric. 

Let 

H = H 2 + H 3 + H 4 + --- + H 2n (2.53) 

Substituting Eqs. (2.41), (2.43), and (2.51) into Eq. (2.50) and rearranging based on 
the order of the state variable, we get 


#2 = 2 ? 2 ( X ) 


1 dV 2 (x) 

2 dx 


BR -1 B 


dV 2 (x) 


dx 


+ 


dV 2 (x) 

dx 


Ax = 0 


(2.54) 



m 




^(x) 


tT 


dx 


1 dV 3 (x) 

2 dx 


BR -1 B t 


5F 2 (x) 


1 T 


dx 


+ a^) Ax+ a^) /2(x) = 0 (2 55) 


#4 = 




9V 4 (x) 


dx 


1 dV 3 (x) 

2 dx 
dV 4 (x) 


BR _1 B t 



T 

!^ 4 ( x ) ! T 

dl 2 (x) 

dx 


! & BE B 

dx 


dx 


Ax+ «!P /!(x) + £}p /3(x) 


dx 


dx 


(2.56) 


*2» = 


l ^x ), 

2 dx 


9y 2n (x) 


-lT 


dx 


1 dV 2n (x) 


2 dx 
dV 2n -i(x) 


BR -1 B r 
/ 2 (x) + b 


dV 2 (x) 


dx 

dV 2 (x) 


+ ^A* 


dx 

/ 2 n-l(x) = 0 


dx — ■ 1 dx 

Solving Eq. (2.54), the algebraic Riccati equation is obtained, 

A T Pi,i + Pi,iA + Q M - P 1)1 BR- 1 B r P ljl = 0 


(2.57) 


(2.58) 


Solving the algebraic Riccati equation, and substituting P^i back into Eq. (2.55), 
P 2) i can be found. The remaining values of the Lyapunov matrix can be found by 
successive substitution of the lower order Lyapunov parameters into higher order Lya- 
punov parameters until P nn is found. After all Lyapunov parameters have been 
determined, stability of the Lyapunov function has to be checked: if F(x) > 0 and 
q{x) > 0, then the nonlinear optimal control is asymptotically stable in the large. 


2.2.3 Nonlinear Optimal Control Procedure 

The general procedure for solving the Lyapunov function V(x) is as follows. 

Step 1: Find the highest order (2n-l) of the state variable from the state equation, so 
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that n can be properly selected. 

Step 2: Assume q(x) and V(x) have the form of Eqs. (2.43) and (2.51) respectively. 
Step 3: Solve the algebraic Riccati equation from Eq. (2.58) and differentiate V^(x) 
and U 3 (x) with respect to x to find dV 7 {x)/dx and dV z (x)/5x. 

Step 4: Substitute dV 2 (x)/dx and dV 3 (x)/dx into Eq. (2.55) and solve it to obtain 
the values of matrix P 2 ,i- 

Step 5: Check the positive definiteness of the Lyapunov function. If it is not stable, 
then adjust the performance index and repeat step 4. 

Step 6: Repeat the process until V 7n (x) is determined. 

The outline procedure will always result in a positive definite V(x). Note that for 
stability V^(x) has to be positive definite. 


Feedback Using Instantaneous Control 

Now' we are in a position to apply the above methodology to our problem and solve 
the equation of motion for linear as w r ell as nonlinear feedback cases. Let us rewrite 
our equation of motion as 

<p + <j) — + [i 7 (j^ <j> + b 7 <f)^ i" u (2.59) 

The additional term which we are adding here is u, the control feedback input, w'hich 
was not used when we discussed open loop system.lt is to be noted that u is scalar in 
this case. 

Defining the state vector x = (x\,x 7 ) T = (<f>, <j>) T , the Eq. (2.59) can be written 
in the state variable form as 

Xi = X 2 (2.60) 

x 2 = -w 2 xi + fjii 2 + /(x) + u . (2.61) 

w'here 



(2.62) 


/(x) = bix\ + n 2 x\x 2 + kxixl 

If Eqs. (2.60) and (2.61) is written in the form as Eq. (2.40), then 


0 1 


Mi 


(2.63) 


B = (2.64) 

1 J 

The method of finding optimal control for this nonlinear model starts by choosing the 
highest order term of / (x) i.e., 2n — 1 = 3 or n = 2. Therefore the performance index 


becomes 


1 f°° I x Ql,l Ql,2 x 


1 r°° 

~ 2 Jo 

i. 


-ft r > dt 


[KJ LQm J 

Note that R=1 here. Now the Lyapunov function is assumed to be 


(2.65) 


V(x) = 


1 X Pj.1 Pi, 2 


Setting 


2 x 2 P 2)1 P 2 , 2 x 


Qi,i — 


and solving the algebraic Riccati equation, we get 


Pi,i = 


P U (1,1) Pi,i(l,2) 
Pi, i (2,1) P u (2,2) 


where 


Pi,i(l, 1) = w 2 Mi + 2^ + ~ + + + 

Pi,i(l, 2) = Pi,i(2, 1) = -V + (u 4 + 1) 1/2 
Pi,i(2, 2) = mi + - 8a; 2 + 8(c 4 + l) 1 ' 2 + 4] 1 ' 2 


( 2 . 66 ) 


(2.67) 


( 2 . 68 ) 


(2.69) 

(2.70) 

(2.71) 


^(x) = ^X^PjjX 


-Ti- 



Therefore, 


rsF 2 (x)i 

[ dx 

T 

= Pl,lX = 

Pi,i(l, 1)^ +P M (1,2 )x 2 



. P M(2,l)ar 1 + P 1>1 (2, 2)x 2 


and 

V>(x) = 5* Tp w* 2 + i ( * J )’> 14 x 

Therefore 

_ r 

3P 2) i( 1, l)xf + 2P 2)1 (2, l) Xl x 2 + P 2 , i(l,2)x| 
P 2 ,i( 2, l)xf + 2P 2 j1 (1 j 2)x iX2 + 3P 2>1 (2, 2)x\ 
Substituting Eqs. (2.72) and (2.73) into Eq. (2.55) and solving we get 

P 2,l = Pi, 2 = 

where 


raH(x)| 

T 

dx J 



Ai(U) p 2 ,i(i,2) 

A, (2,1) P 2tl ( 2,2) 


A,i(l,2)=P 2J (2,l) = feoO; !) 

I“ 2 + Pi,i(I,2)J 


(2.72) 


(2.73) 


(2.74) 


(2.75) 


P2>1 (1 ’ 1} 3 [ Q2 ' 1 (1 » 2 ) + P m( 2 » 2)P 2i i(1, 2) + 2P 11 (1 ) 2)P 2i 1 (1, 2) + 2u; 2 P 2)1 (l, 2) 


-^iA,i(1,2)] 


(2.76) 


ft.i(2,2) = 2) - P 2-I (i,2)- A (i, 211 , 

3^+3^2) (2-77) 

Since n = 2, we have 


'«(*) = ^(x 2 ) t P 2 , 2 x 2 


and 


pF 4 (x)‘ 

T 

dx 



2^iP 2 , 2 (1, 1) + 2xix 2 P 22 (1 ) 2) 

2^ 2 P 2)2 (2, 2) + 2x 2 xfP 2>2 (l, 2) J 


substituting Eqs. (2.72), (2.73) and (2.78) into Eq. (2.57), we get 


P 2,2 = 


J\ 2 (l,l) P 2 , 2(1,2) 
^,2(1,2) P 2)2 (2, 2) 


(2.78) 


(2.79) 


P 2.2(1, 2) = P 22 (2. 1) 





„ „ _ 2[P 2 , 1 (l,2)] 2 + 3P ! , 1 (2,2)P 2il (l,2)-(3 M (l,2)-6 2 P 1 ,,(l,2)-,/2fi,,(2,2) 

Pw(1 ’ 2) ^i-2P;,(2,2) 

(2.80) 

P2,2<1, 1) = [P 2 ,i(l, 2)] 2 + w 2 P2,j( 1, 2) - 0.5ftP M (l,2) + P,.,(l, 2)P 2 , 2 (1, 2) (2.81) 

P [6P2,i(2,2)P2.,(l,2)-2P 2i (l,2)-6 1 P lil (l,2)-62Pi,i(2,2)l 

P22<2 ' 2) 2[oj 2 + Pi.i(l,2)] (282) 

On substituting Eq. (2.72) into Eq. (2.49) we get the linear feedback optimal as 


u(x) = -[Pi,i(l, 2)x x + Pi,i(2,2)x 2 ] (2.83) 

Similarly on substituting Eqs. (2.72), (2.73) and (2.78) into Eq. (2.49) we get the 
nonlinear feedback optimal control as 


u(x) 


-[Pi,i(l, 2)x x + P M (: 2, 2)x 2 + 3xlP 2jl (2, 2) + + 2x 1 x 2 )P 2)1 (l, 2) 

+2x|P 2)2 (2, 2) + 2x 2 x\P^(l, 2)] (2-84) 


Furthur, these equations of u(x) have been substituted in the equation of motion and 
Eqs. (2.60) and (2.61) have been solved with the help of Runge-Kutta method of 
4 th /5 th order on MATLAB to obtain plots of roll angle and roll rate against time. 


Feedback Using an Actuator for Control Surface 


In the above case, we have assumed an imaginary actuator which is applying the 
feedback instantaneously as the output is obtained. In practical situation, such an 
assumption cannot be made. Therefore, in this part a modification in the above 
problem has been made by adding one more differential equation which will take care 
of the above discrepancy. 

An aileron control surface with first order actuator dynamics is modeled as 

&a — —~Sa + -&Acom - (2.85) 

r r 

where 6a is the aileron angle, 6a C om * s commanded aileron angle and r is the time 
constant of the actuator, i.e.. feedback time lag. 



Now defining our state vector as x = (x 1) x 2 ,x 3 ) r = (4,i>,5 A ) T , our equation of 
motion can be written in the state variable form as 


Xl = x 2 


( 2 . 86 ) 

(2.87) 

( 2 . 88 ) 


X 2 --UIJ + Ml z 2 + /( x ) + Xi 
x 3 = -kx 3 + ku 

where k = K If Eqs. (2.86), (2.87) and (2.88) are written in the form as Eq. (2.40), 
then 

0 1 0 
-w 2 Hi 1 

0 0 -k 


A = 


(2.89) 


and 


B = 


0 

0 

k 


(2.90) 


Setting 


1 0 0 
Qi.i = o l o 
0 0 1 

and solving the algebraic Riccati equation we get 


(2.91) 


Pi, i(M) PiA 1,2) Pi,i(l,3) 

Pi,i = I Pu (2,1) Pi,i(2,2) P u (2,3) (2-92) 

Pi,i(3, 1) Pi,i(3,2) Pi,i(3, 3) 

As in the previous case, we can substitute the expression for <?V2(x)/9x. dV 3 (x)/dx 
and 

dV 4 (x)/cbc in Eqs. (2.55) and (2.57) and solve these equations to obtain P 2.1 and P 2j2 
where 

"Pa, i(l,l) P 2)1 (l,2) P 2 ,i(1,3) 

P 2) i=Pi, 2= | P 2ll (2,l) P 2 ) i(2, 2) P 2 ,i(2, 3) | (2.93) 

P 2 .i(3.1) P 2 ,i ( 3,2) P 2> i (3, 3) 



and 


^ 2 , 2 ( 1 , 1) P 2>2 (1,2) P 2 , 2 (l,3) 

P 2 ,2= | P 2>2 (2, 1) P 2j2 (2, 2) P 2 , 2 (2,3) | (2-94) 

P 2 , 2 (3,1) P 2 , 2 (3,2) P 2 , 2 (3,3) 

Although all the elements of the matrix P 2> i and P 2 2 can be obtained by solving 
Eqs. (2.55) and (2.57), only those elements of P 2) i and P 2j2 are given below which 
are required to solve our state variable equation. 

Qi,2(1. 3) + [Qi < 2 (1, l)/a> 2 ] 


Pi.2(1,3) = 


2FP lil (3, 1) + k 2 Pi t i(3 , 3 ) + k + [A: 2 P li i(3, l)/u; 2 ] 
@i, 2 (l. 1) _ fc 2 Pi,i(l> 3 )Pi )2 (1, 3) 


Pi, 2 (1,2) = 




P 1 , 2 (3,3) = 


n .0 _ Pb(l,3)Il-*»A4(3,2)] + ^(l,2) 

Pl ' 2(2,3) ^ + ««u(l,3) 

@ 1 , 2 ( 1 , 3) - 2fc 2 P 1 , 1 (3,3)P 1>2 (l,3) -^,2(2,3) - 2fcP 1>2 (l,3) 

3fc 2 P u (3,l) 

r, f ... [Pi, 2 (1,3)3 2 

P2 ’ 2(1 ’ 3) "-^iT 

p &2Pi,i(3,3)-2fc 2 P 1 ,2(2,3)P 1 , 2 (l,3) 

2fc 2 P u (3, 1) 


P 2 , 2(3,3) = - 


3Pi, 2 (1.3)P 1)2 (3,3) 


Pi, i (3,1) 

The linear feedback optimal control for this case is 

u(x) = -fc[P u (l,3)xi + Pi,i (2, 3)x 2 + Pi.i(3,3)x 3 ] 


(2.95) 

(2.96) 

(2.97) 

(2.98) 

(2.99) 
( 2 . 100 ) 
( 2 . 101 ) 

(2.102) 


and the nonlinear feedback is 


u(x) = — fc[P 1)1 ( 1, 3)xi + Px i i(2, 3)x 2 + Pi,i(3,3)x 3 + Pi, 2 (1,3)x 2 + Pi, 2 (2,3)x 2 
+2Pi, 2 (1,3)xxX 3 + 2Px >2 (2, 3)x 2 x 3 + 3Px ) 2(3, 3)x 3 
+2P 2>2 (1,3)x 2 x 3 + 2P 2)2 (2,3)x 2 x 3 4- 2P 2 , 2 (3, 3)x^] (2.103) 

Substituting the values of u(x) in Eq. (2.88) and solving Eqs. (2.86), (2.87) and (2.88) 
with the help of Runge-Kutta method of 4 th /5 th order on MATLAB, we obtain the 
plots of roll angle, roll rate and aileron angle against time. 



Feedback Considering Sideslip Effects 

Until now, we have not considered the effects of sideslip on the motion of aircraft. In 
this section we will see the variation of sideslip angle and sideslip rate with time and 
their overall effect on other parameters. 


We define our state variable vector as x = (a*, x 2 , x 3 , x 4 , x 5 ) T = (<£, <5 A , /?, j3) T . 


Our equation of motion can be written as 

x'i = X2 (2.104) 

x 2 = —io 2 X\ + H\X 2 + L$x a - L r x$ + /(x) 4- x 3 (2.105) 

x‘ 3 = — fcx 3 + ku (2.106) 

x 4 = X5 (2.107) 

x 5 = —NpX 4 + N r x 5 - N p x 2 (2.108) 


where, L@ is the dihedral effect, L r is the roll moment due to yaw rate, N$ is the 
directional stability, N r is the yaw rate damping and N p is the yaw moment due to 
roll rate. 

For this case the corresponding A and B matrices are 



1 0 

Mi 1 
0 -k 
0 0 
N p 0 


0 0 

Lp —L r 

0 0 

0 1 

-N 0 N r 


B = 


0 

0 

k 

0 

0 


(2.109) 


( 2 . 110 ) 
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Setting 


1 0 0 0 0 

0 10 0 0 

Qm= 0 0 1 0 0 (2.111) 

0 0 0 1 0 

0 0 0 0 1 

Using these values of A, B and Qi,i and solving the algebraic Riccati equation we get 
Pi,i where 

P h i(l.l) A,i(1,2) P i,i (1,3) P u ( 1,4) P u (l,5) 
p h i(2,l) A, i(2,2) Pi,i(2, 3) P M (2,4) P 1(1 ( 2,5) 

Pi,i(3,l) Pi,i (3, 2) P M ( 3,3) Pi,i(3, 4) P M (3,5) (2.112) 

Pu(4,l) Pi,i(4, 2) Pi, 1 (4,3) P u (4,4) P lil (4,5) 

Pi,i(5, 1) Pi,i (5, 2) P 1,1 (5, 3) Pi,i(5,4) P a)1 (5,5) 

Now solving the Eqs. (2.55) and (2.57) after substituting the expressions of dV^x) /dx. 
dV 3 (x)/dx and dV 4 (x)/dx in them we get P 2) i and P 2)2 where 


P M = 


P 2,1 = P 1.2 = 


1,2 


P 2)1 (l,l) P 2 ) i(1,2) P 2 ,i( 1, 3) P 2>1 (1,4) P 2)1 (l,5) 

P2.i(2,1) P 2) i(2, 2) P 2|1 (2, 3) P 2>1 (2,4) P 2>1 (2,5) 

P 2 ,i(3,1) P 2i i(3,2) P 2) i (3, 3) P 2>1 (3,4) P 2iI (3,5) 

Pa, i (4,1) P 2 > i(4,2) P 2 ,i(4, 3) P 2i1 (4,4) P 2 ,i(4,5) 

P2,i( 5,1) P 2 ,i(5, 2) P 2 ,i(5, 3) P 2 ,,(5,4) P 2i1 (5,5) 


(2.113) 


and 


P 2l2 (l,l) P 2 , 2 (l, 2) P 2 , 2 (l,3) P 2 j2 (1,4) P 2>2 (1, 5) 

P 2>2 (2, 1) P 2 , 2 (2, 2) P 2 , 2 (2,3) P 2i2 (2,4) P 2 , 2 ( 2,5) 

P 2 , 2 = | P 2 , 2(3,1) P 2)2 (3, 2) P 2)2 (3,3) P 2)2 (3, 4) P 2 , 2 (3,5) | (2.114) 

P 2)2 (4, 1) P 2 , 2 (4, 2) P 2j2 (4, 3) P 2 , 2 (4,4) P 2 , 2 (4,5) 

P 2)2 (5,l) P 2 j2 (5, 2) P 2>2 (5,3) P 2 , 2 (5,4) P 2j2 (5,5) 

Those values of matrices P 2;1 and P 2)2 are given below which required to solve our 
equation of motion. 

^ 2 <3i. 2 (1, 3) + Qi, 2 (1, 1) 


Pi. 2 (3,1) = 


w 2 [2FP u (3, 1) + 2A' 2 Pj il (3.3) + k} + fc 2 P u (3, 1) 


(2.115) 
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P 1,2(3, 3) 


A. 2 ( 3,4) = 

d „ «_'Pi.i(3,5)P u (3,4) 

PU(3 ’ 5) Pu(3,4) 

Pi j(2j 1} = Qu(l,l)-^Pu(3,I)P u (3.1) 

’ ’ a; 2 

P "> ox - [2-2fe 2 P 1 , 1 (3,2)]P 1 , 2 (3,l) + 2F 1 , 2 (2,l) 
h2{ ’ ’ 2A: 2 P 1) j( 3, 1) + 2w 2 

Qi, 2 (3, 1) — 2P 1)2 (3, l^P^, 3) + fc] - a; 2 P x , 2 (3, 2) 


P 2 , 2 (3,l) 


3fc 2 Pi,i(3,l) 

[Pi, 2(3, l)] 2 

Pi, i (3,1) 


p ^ &2Pi,i(3,3) - 2fc 2 P 1 , 2 (3,2)P 1)2 (3,l) 

2 ’ 2 ^ 3, / 2A: 2 P 1 , 1 (3,1) 


P 2 , 2 (3,3) = - 
P 2>2 (3, 4) = 
P 2 , 2 (3, 5) = 

The linear feedback in this case is 


Pi, i (3,1) 

Pi, 2 (3, 1)Pi, 2 (3,4) 
Pi, i (3,1) 

Pi, 2 (3, 1)Pi, 2 (3, 5) 
Pi, i (3,1) 


(2.116) 

(2.117) 

(2.118) 

(2.119) 

( 2 . 120 ) 
( 2 . 121 ) 
( 2 . 122 ) 

(2.123) 

(2.124) 

(2.125) 


tt(x) = — fc[Pi,i(3, l)xi+Pi,i(3, 2 )x 2 +Pi > i(3,3)x 3 +Pi i i(3,4)x 4 + P i,i(3,5)x 5 ] (2.126) 


and the nonlinear feedback is 

u(x) = -^1,1(3,!)^ +P 1 ,i(3,2)x 2 + Pi !l (3,3)x 3 + P M (3,4)x 4 + P!,i(3,5)x 5 

+Pi, 2 (3, 1)x 2 + P 1>2 (3, 2)x 2 + 2P 1i 2 (3, l)xiX 3 + 2P li2 (3, 2)x 2 x 3 + 3P 1 , 2 (3 > 3)x 2 
+2 Pi, 2 (3, 4)x 3 x 4 + 2P 1i2 (3, 5)x 3 x 5 + P u (3, 4 )x^ + Pi, 2 (3, 5)x! + 2P 2 , 2 (3, l)x 2 x 3 
+2P 2 i2 (3, 2)x\x 3 + 2P 2>2 (3, 3)x^ + 2P 2>2 (3, 4)x 3 x^ + 2P 2>2 (3, 5)x 3 x^] (2.127) 

Finally these values of u(x) are substituted in the state variable equations and solved 
with the help of Runge-Kutta method of 4 t/l /5 th order on MATLAB to obtain plots 
of roll angle, roll rate, aileron angle, sideslip angle and sideslip rate against time. 



Chapter 3 


NUMERICAL RESULTS 


For the demonstration of the methodology discussed in the previous Chapter, same 
delta wing, which has been used in Refs. [2] and [8], is taken here. The particulars of 
the delta wing is given below. 

Leading-edge-sweep angle = 80 deg 

Aspect ratio = 0.705 

Chord = 0.429m 

L c = Chord/4 

Area = 0.0324m 2 

Mass = 0.284kg 

I x = 0.27 x 10~ 3 kg-m 2 

Pair = 0.12 x 10~ 2 kg/m 3 

U c = 15m/s 

p x = 0.378 x 10~ 4 kg-m 4 

This value of p x was found by trial and error to be the only value for which onset angle, 
period and amplitude agree with the experiment. These choices lead to C\ = 0.354 
and C 2 = 0.001. 
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3.1 Response of Open Loop System 

To know about the character of motion at non-equilibrium positions, we solve Eqs. 
(2.28) and (2.29) using Runge-Kutta method of 4 th /5 th order on MATLAB. The plots 
which have been obtained, of roll angle and roll rate with respect to time, indicates 
the wing rock motion to be a limit cycle motion. 

The amplitude of the limit cycle motion depends on the angle of attack as it can 
be observed from Fig. (3.1). An early transition to limit cycle motion along with 
increased amplitude is observed, as the angle of attack increases. The corresponding 
phase plane plots for second order system of roll angle and roll rate for three angles 
of attack is shown in Fig. (3.2). Fig. (3.3) show’s the w'ing rock motion triggered by 
roll angle for fifth order system at a = 25 deg. Fig. (3.3a) and Fig. (3.3b) exhibit 
limit cycle motion for roll angle and roll rate respectively whereas Fig. (3.3c) and 
Fig. (3.3d) exhibit limit cycle limit cycle motion for sideslip angle and sideslip rate 
respectively. It can be observed from Fig. (3.3d) that initially sideslip rate increases 
sharply and then transits into limit cycle motion. Fig. (3.3e) show’s the corresponding 
phase plane plot of roll angle and roll rate whereas, Fig. (3.3f) show’s the phase plane 
plot of sideslip angle and sideslip rate. 

It has been stated earlier that limit cycle motion may be initiated by sideslip as 
well. In fact, this behaviour has been observed, even for very small sideslip angle of 
0.01 rad. Fig. (3.4) show's the wing rock motion triggered by sideslip angle for fifth 
order system at a — 25 deg. Fig. (3.4a) and Fig. (3.4b) show the roll angle and roll 
rate, starting from zero, attain limit cycle motion in 1200 sec. Fig. (3.4c) shows that 
the sideslip angle, starting from 0.01 rad, after initial fluctuations attains limit cycle 
motion. Similar behaviour has been observed from Fig. (3.4d) for sideslip rate, but in 
this case its initial value is 0 rad/s. Fig. (3.4e) shows the corresponding phase plane 
plot of roll angle and roll rate and Fig. (3.4e) and Fig. (3.4f) shows the corresponding 
phase plane plot of sideslip angle and sideslip rate. This exhibition of wing rock by 
sideslip is due to the small dihedral effect which has been considered in our model. 
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3.2 Response of Closed Loop System 


For the second order state variable, we observe that our model gets stabilized in 
approximately 6 secs, for small initial conditions of (<f>,4>)={ 0.35 rad, 0 rad/s) for 
both nonlinear and linear feedback cases. In fact, no appreciable difference has been 
observed between nonlinear and linear feedback responses for these initial conditions 
as it can be observed from Fig. (3.5). Fig. (3.5a) shows nonlinear feedback response 
at a = 25 deg whereas Fig. (3.5b) shows linear feedback response at o = 25 deg. 
Fig. (3.5c) shows the corresponding phase plane plot of roll angle and roll rate of 
nonlinear and linear feedback response. In Fig. (3.3a) and Fig. (3.3c) the results 
obtained are compared with the results obtained in Ref. (14). As it can be observed 
from these plots that the results obtained in this thesis is same as what it has been 
obtained in Ref. (14). Also it has been observed for nonlinear as w T ell as linear cases 
that for these small initial conditions feedback responses at different angles of attack 
are found to be same. Whereas for large initial conditions of (<£,<£)=(1.4 rad, 3.5 
rad/s), nonlinear feedback response proves to be better than linear feedback response, 
as it can be observed from Fig. (3.6). For a = 25 deg, linear feedback response even 
gets diverged in 0.37 seconds. But in case of wing rock, the amplitudes of limit cycle 
motion for our model at a = 25 deg are (<^>, <^>)=(0.6 rad, 0.08 rad/s). Thus at these 
initial conditions, it has been observed that both nonlinear as well as linear feedback 
responses are found to be equally effective in supressing wing rock motion. Here it 
has been assumed that 1) = 1/64, Qi^i 1, 2) = 0 and < 22 , 2 ( 1 , 2 ) = 0. 

For the third order state variable, similar features like second order have been 
observed. But in this case, for initial conditions of (4>,^8a) =(0.35 rad, 0 rad/s, 0 
rad), the time to stabilize has increased to 8 secs., as it can be observed from Fig. 
(3.7a) and Fig. (3.7b). This expected increase in time is due to the use of actuator 
with time constant r = 0.0495 s. The maximum aileron deflection has been observed 
to be -0.21 rad for both nonlinear as well as linear feedback response. Fig. (3.7c) 
show r s the corresponding phase plane plot of roll angle and roll rate. \\ hen the initial 
conditions are {<t>,4>,8 A )=( 0.6 rad. 0.08 rad/s. 0 rad), the stabilization time remains 



almost same, but the maximum aileron deflection has increased to -0.45 rad. as it can 
be observed from Fig. (3.8a). Fig. (3.8b) shows the phase plane plot for nonlinear 
and linear feedback response at a = 25 deg. and Fig. (3.8c) shows the comparison 
phase plane plot of nonlinear feedback response at different angles of attack. It has 
been assumed that (? u (l, 1) = 1/64 and <?„( 1, 3) = 0. The P u matrix, which has 

been obtained by solving the algebraic Riccati equation for three angles of attack are 
given below. 

1.97 1.4404 0.048 

1.4404 2.8934 0.098275 

0.048 0.098275 0.02386 

1.9734 1.4476 0.0483 
1.4476 2.8889 0.0981 
0.0483 0.0981 0.0239 

1.9751 1.451 0.0484 

1.451 2.8875 0.0981 

0.0484 0.0981 0.0239 

For the fifth order state variable, where we have also considered beta dynamics, the 
nonlinear and linear feedback responses for the initial conditions (</>, </>, 5 A , /?, /?)=(0.35 
rad, 0 rad/s, 0 rad, Orad, 0 rad/s) is found to be almost overlapping on each other. 
But for the clarity, only nonlinear feedback response is shown in Fig. (3.9). In Fig. 
(3.9a) responses of all the five state variables have been shown with respect to time. 
Fig. (3.9b) shows the corresponding phase plane plot of roll angle and roll rate at 
different angles of attack and Fig. (3.9c) shows the corresponding phase plane plot 
of sideslip angle and sideslip rate at different angles of attack. Fig. (3.10) show’s the 
nonlinear and linear feedback response for the initial conditions (<£,<£,5x, /?,/?)= (0.6 
rad, 0.08 rad/s, 0 rad, 0 rad, 0 rad/s). In this case the roll angle, roll rate and aileron 
angle stabilize in almost 8 sec, but the sideslip angle and the sideslip rate continue 
to fluctuate till 45 sec. Fig. (3.11) shows the corresponding phase plane plots. Fig. 
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(3.11a) shows the phase plane plot of roll angle and roll rate for nonlinear and linear 
feedback response at a = 25 deg. whereas Fig. (3.11b) shows the corresponding phase 
plane plot for sideslip angle and sideslip rate. A comparison phase plane plot for 
nonlinear feedback response at different angles of attack for roll angle and roll rate 
and sideslip angle and sideslip rate has been shown in Fig. (3.11c) and Fig. (3.11d) 
respectively. Very little variation in these phase plane plots have been observed as the 
angle of attack varies. It has been assumed that the magnitude of L 0 = -0.02822, 
L r = 0.1517, Np = 1.3214, N r = —0.2491 and N p = —0.0629, for all the three angles 
of attack considered in this thesis, because of their negligible variation for these three 
angles of attack. In this case also Q 1 j 2 (1, 1) = 1/64 and Q 1 j 2 (1,3) = 0 has been 


assumed. The P^i matrix for three angles of attack are given below. 


Px 


>*0=25 deg 


Pi4a=: 


22.5 deg 


1 Ma=21.5 deg 



1.9726 

1.4408 

0.048 

-0.0726 

-0.1504 


1.4408 

2.9118 

0.0989 

0.2209 

-0.372 

= 

0.048 

0.0989 

0.0239 

0.0078 

-0.0126 


-0.0726 

0.2209 

0.0078 

4.7865 

0.3642 


-0.1504 

-0.372 

-0.0126 

0.3642 

3.5666 


1.9759 

1.448 

0.0483 

-0.0709 

-0.1506 


1.448 

2.9067 

0.0987 

0.221 

-0.3711 

= 

0.0483 

0.0987 

0.0239 

0.0078 

-0.0125 


-0.0709 

0.221 

0.0078 

4.7867 

0.3642 


-0.1506 

-0.3711 

-0.0125 

0.3642 

3.5667 


1.9776 

1.4514 

0.0484 

-0.0702 

-0.1507 


1.4514 

2.9053 

0.0987 

0.2210 

-0.3707 

= 

0.0484 

0.0987 

0.0239 

0.0078 

-0.0125 


-0.0702 

0.221 

0.0078 

4.7867 

0.3742 


-0.1507 

-0.3707 

-0.0125 

0.3642 

3.5667 


A comparison between second, third and fifth order state variable has been shown 
in Fig. (3.12). From this figure, it can be observed that the stabilization time for 
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the third and fifth order state variable remains approximately same. This is due to 
the fact that same actuator has been used for the two cases and the dihedral effect 
considered here is quite small. If dihedral effect is large, then appreciable difference 
between the third order and the fifth order state variable will be observed, even if 
same actuator is used for both the cases. 


Fig.(3.1a) Limit cycle motion of Roll angle at alpha=21 JS deg. 


Timet*). 

Rg.(3. lb) Limit cycle motion of Roll tale at ahsha=2lJ5 deg. 


Fig.(3.1c) Limit cycle motion of Roll angle at alpha=22.5 deg. 


Hg.(3.1d) Limit cycle motion of Roll rate at alphass22.5 deg. 


Fig.(3.1e) Limit cycle motion of Roll angle at aJpha=25 deg. 


Fig.(3Jf) l imit cycle motion of Roll *mc at aJ$***25 deg. 


Figure 3.1: Comparison of limit cycle motion triggered by roll angle it different angles 


of attack for second order state variable. 





Roll angle (rad) 

Fig.(3.2c) Phase plane plot for open loop system at alpha=25 deg. 


Figure 3.2: Comparison of phase plane plots at different angles of attack for second 
order state variable. 
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Fig.(3.3e) Phase plane plot for open loop system at alpha=25 deg. 





Rg.(3.30 Phase plane plot for open loop system at alpba=25 deg 


Figure 3.3: Wing Rock motion triggered by roll angle for fifth order state variable at 
a = 25 deg. 
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Roll rate (rad/s) Sideslip angle (rad) ang,c 




Hg.(3.4c) Limit cycle motion of Sideslip angle at alpha=25 deg. 



- 0.2 0 02 
Roll angle (rad) 

Fig.(3.4c) Phase plane plot for open loop system at aJpha=25deg. 


Figure 3.4: Wing Rock motion triggered by j 
at a = 25 deg. 
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Sideslip angle (nd> 

Fig.(3.4f) Phase plane pkx for open loop system at alpha-25deg. 


ip angle for fifth order state variable 
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Fig.{3.5a) Nonlinear feedback response at alpha=25 deg. 
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Figure 3.5: Nonlinear and linear feedback response for second order state variable at 
a = 25 deg, with initial conditions (<?>, 0)=(O.35 rad, 0 rad/s). 
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Figure 3.6: Comparison of nonlinear and linear feedback response for second order 
state variable at different angles of attack, with initial conditions (<fi, <j>)={\A rad. 3.5 
rad/s). 
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Figure 3.7: Nonlinear and linear feedback response at a — 25 deg for third order state 
variable with initial conditions rad, 0 rad/s, 0 rad). 
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Figure 3.9: Nonlinear feedback response for fifth order state variable with initial con- 
ditions (<f>,j>,6 A ,(3,0)=( O.35 rad, 0 rad/s, 0 rad, 0 rad, 0 rad/s). 


46 





Sideslip angle (rad) , Roll *igle (rad) 







Figure 3.10: Nonlinear and linear feedback response for fifth order state variable at 
a = 25 deg with initial conditions (<j>, <f>, S A , /?, 0) =(0.6 rad, 0.08 rad/s, 0 rad, 0 rad, 
0 rad/s). 
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considering sideslip effect at alpha=25deg 



considering sideslip effect at alphas25deg 




sideslip effect at different angles of attack sideslip effect at different angles of attack 


Figure 3.11: Phase plane plots of closed loop system for fifth order state variable with 
initial conditions ( <f > , /?,/?)= (0.6 rad, 0.08 rad/s, 0 rad, 0 rad, 0 rad/s). 
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Fig.(3.12b) Nonlinear feedback response for Roll rate at alpha=25deg 



Figure 3.12: Comparison of three variables at a = 25 deg as their order of state 
changes. 
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Chapter 4 

CONCLUSIONS 

A scheme has been applied for controlling the wing rock motion of aircraft flying at 
angles of attack near stalling onset, where optimal solution of nonlinear system is 
obtained by using Lyapunov function and. performance index, which are both positive 
definite. Earlier, this method had been used for obtaining the solution of second order 
state variable system, but in this thesis, it has been extended to third and fifth order 
state variable system. Both open loop and closed loop systems have been considered. 
The plots obtained for open loop system have demonstrated limit cycle motion of the 
model when it is not in equilibrium position. The phase plane plots have shown that 
the limit cycle amplitude depends on angle of attack and not on initial conditions. It 
has been observed that sideslip angle can also initiate wing rock motion if dihedral 
effect considered is quite small. For the closed loop system, both nonlinear and linear 
feedback responses have been obtained. From these plots, we can observe that both 
nonlinear and linear feedback controls are equally effective for supressing the wing rock 
motion of our model. According to the expectation, the actuator considered in third 
and fifth order state variable system increases the stabilization time of the model. 
As observed for the fifth order state variable system, roll angle, roll rate and aileron 
angle stabilize in very less time, but the disturbance in sideslip angle and sideslip 
rate, although small in magnitude, stays for quite long time. It is to be noted that 
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in Ref. [12], it has been assumed that (j> — 2/3, but from the plots obtained for our 
model, no such relation has been observed. This difference can either be due to wrong 
assumption or due to different model considered in Ref. [12]. Whatever ever be the 
reason, but it is certain that such assumption cannot be considered blindly. 

In this work, variations in flight conditions have not been considered. So, consider- 
ing the practical aspect of the work, if there is a provision made which can update the 
coefficients of the aircraft as the flight condition changes, then these optimal control 
equations can be utilised for controlling the wing rock motion of the actual aircraft. 
Future scope of the work is that, the nonlinearities in the actuator can also be con- 
sidered which have been neglected in this thesis. It has been assumed that the wings 
remain straight during the flight, which may not be so especially when the aircraft is 
undergoing wing rock motion. So, the structural stiffness of the wings may also be 
taken into account, which has been neglected in this work. 
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